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Description 



Background of the Invention 
Field of the Invention 

The present invention relates to velocity measurement systems and, more particularly, to such systems applying 
correlation methods. 



Description of the Related Technology 



There is a class of systems that measure distance or velocity with sound waves. These systems include underwater 
sonar and medical imaging equipment (hereinafter generally referred to as "sonar systems"). Underwater sonar sys- 
tems, in particular, may be manifested as current profilers, velocity logs, and so forth. 

A current profiler is a type of sonar system that is used to remotely measure water velocity over varying ranges. 
Current profilers are used in freshwater environments such as rivers, lakes and estuaries, as well as in saltwater 
environments such as the ocean, for studying the effects of current velocities. The measurement of accurate current 
velocities is important in such diverse fields as weather prediction, biological studies of nutrients, environmental studies 
of sewage dispersion, and commercial exploration for natural resources, including oil. 

Typically, current profilers are used to measure current velocities in a vertical column of water for each depth "cell* 
ot water up to a maximum range, thus producing a "profile" of water velocities. The general profiler system includes a 
transducer to generate pulses of sound (which when downconverted to human hearing frequencies sound like "pings") 
that backscatter as echoes from plankton, small particles, and small-scale tnhomogeneities in the water. 

Two important types of sonar technology that may be used for current profiling are Doppler and correlation. The 
basic concept of velocity measurement using signal correlation has been known for many years. For example, in 1 964, 
Dickey, Jr. was issued U.S. Patent No. 3,147,477 ("Dickey '477") disclosing a system to measure speed including a 
source of wave energy (or transmitted signal) directed toward a body, multiple receiving devices having a known sep- 
aration for receiving echoes from the body, correlating means including time delay means, and speed determining 
means dependent on the quotient of the receiver separation and the time delay. Although the patent document describes 
an embodiment of a correlation system for an aircraft using radar (electromagnetic signals), it suggests a correlation 
system for a water-going vessel using sonar (sound signals). 

A correlation system inherently requires that the transmitted signal must have some similarity with itself to achieve 
a maximum correlation at the specified time delay. A correlation system was disclosed by Roeder, et aL (U.S. Patent 
No. 4,103,302) wherein signal similarity was achieved by generating a pulse train, or series of repeating pulses. Ac- 
cording to the disclosure, the use of repeating pulses provides advantages over the continuous wave (CW) system 
described in Dickey '477 including: a simplified transmission source, insensitivity to variations in attitude of the con- 
veying vehicle, and improved accuracy and reliability derived from digital electronics. 

U.S. Patent No. 4,244,026 to Dickey, Jr., ("Dickey '026") corresponding to EP-A-0 010 974 discusses a correlation 
system using a pulse train from which complex correlation values are employed in velocity measurement. According 
to Dickey '026, the method of signal processing the complex correlation values estimates components of location by 
specifically requiring the steps of processing the correlation values as samples of a continuous function of position and 
curve fitting the samples to the magnitude or real part of the continuous, spacial correlation function to provide an 
estimate of the peak value of correlation magnitude. Thereafter, the location components at the correlation peak are 
divided by twice the pulse repetition period of the pulse train to obtain velocity measurements. This signal processing 
is described in the specification of Dickey '026 at col. 13, lines 7-40. 

As previously noted, correlation sonar technology can be used in current profilers. In designing a current profiler, 
in particular, trade-offs are made among a variety of factors, including maximum profiling range and temporal, spacial 
(the size of the depth cell), and velocity resolution. Temporal resolution refers to the time required to achieve a velocity 
estimate with the required degree of accuracy. In typical applications, a current profiler will make a series of measure- 
ments which are then averaged together to produce a single velocity estimate with an acceptable level of velocity 
variance, or squared error. 

Over the past ten years, acoustic Doppler technology has been successfully developed for short- and medium- 
range current profiling applications such as measuring sediment movement in rivers and estuaries. To achieve the 
longer ranges necessary to study water velocities at deep ocean depths, lower frequencies must be used so that energy 
losses due to sound absorption by the water can be kept at manageable levels. Although there is no fundamental 
limitation in the Doppler technology that prevents the use of lower frequencies, there is a practical constraint due to 
the large size of the acoustic transducers required for these lower frequencies. Existing correlation sonar techniques 
measure current velocities by transmitting two short pulses from a relatively smaller transducer along a correspondingly 
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wider angle acoustic beam. By careful cross-correlation of spatially sampled returns it is possible to determine the 
displacement of the acoustic targets within the beam from one pulse to the next, and hence, compute their average 
velocity. The primary advantage of a correlation sonar is that for comparable sized systems, the correlation sonar is 
able to utilize a lower frequency regime and, therefore, potentially gains substantial current profiling range. 

Acoustic bedload velocity estimates using a broadband pulse-pulse time correlation technique" by D. Sutton and 
J. Jaffe, 1992 ("Sutton") discusses a broadband correlation sonar system using multiple transducers and pulses. In 
Sutton, one receiver time series is autocorrelated in time (i.e., one pulse is autocorrelated to a later pulse). This auto- 
correlation is used in a maximum likelihood algorithm which is capable of finding only one vector component of the 
three-axis velocity vector. This is a fundamental principle of traditional Doppler sonars, in general. In order for the full 
three-axis velocity vector to be determined: it is necessary to have separate beams of sound directed toward (in general) 
different sections of the seabed producing different stochastic time series for each beam. The ML method is then 
applied to each time series in isolation to produce one vector component velocity from each time series in isolation. 
To form the full three -component velocity vector, individual ML estimates are later combined by virtue of the fact each 
one represents a component of velocity in a different direction. 

Nonetheless, it is only recently that the computational burden required of sophisticated signal processing methods, 
required for accurate correlation processing, has been made feasible by advances in the semiconductor and computer 
architecture fields. The improved signal processing should include a means for bottom tracking so that the correlation 
sonar can also serve as a navigation device. It would also be desirable to allow individual users the flexibility of se- 
lectably trading off standard deviation and profiling resolution so as to maximize range. 

Summary of the Invention 

The above-mentioned needs are satisfied by the present invention which includes a correlation sonar system 
utilizing a monopulse transmission method and a maximum likelihood method of signal processing as defined in claim 
1 , 9 and 1 9 : respectively. 

In the presently proposed embodiment of the invention, the correlation system performs the functions of a current 
profiler, which measures current velocities in a water column, and a bottom tracker, which measures vessel velocity 
relative to the bottom of a body of water, such as the ocean. Basically, the present inventive system begins making 
velocity measurements by transmitting a single, uninterrupted acoustic pulse into the water. This pulse is of a long 
duration so as to not only provide a correlation delay that is less than the pulse width, but also so as to typically 
simultaneously ensonify all of the "scatterers" on the bottom. 

The single pulse is typically transmitted for a sufficient duration to encompass both the range of temporal delays 
to all scatterers in the water column cell (and the bottom) as well as the temporal delay for correlation processing. The 
single pulse may be coded, using phase encoding, for example, for the purpose of reducing the noise variance on the 
measured velocity. 

The pulse echoes produced when the pulse encounters small-scale in homogeneities in the water, and then the 
bottom, are received by a set of transducers. The pulse echoes received by the transducers are represented internally 
to the system as complex values. These complex values derived from each transducer are cross-correlated with those 
from other transducers at different time lags (typically zero and tau) in all possible permutations to form complex cor- 
relation values. In this manner, both spatial and temporal information is derived from the process. 

The complex correlation values are not converted to correlation magnitudes but rather the complex correlation 
values, which retain explicit phase information, are used to form a probability function. Hypothesized parameters, such 
as signal-to-noise ratio, seabed slope angles, parameters characterizing seabed roughness, and three component 
velocity (v x , v y , v z ), are combined with the complex correlation values to obtain a probability. The combining of hypoth- 
esized parameters and complex correlation values into a probability is accomplished mathematically using a conditional 
probability density function (for example, see equation B.15 on p. 348 of Statistical Theory of Signal Detection, Carl 
Helstrom, Pergamon Press, 1 968). Hypothesized parameters are then adjusted by a model to maximize the probability 
This method is appropriately known as the maximum likelihood method. 

Typically, the complex values are correlated at zero lag and at the autocorrelation lag where there is a peak in the 
autocorrelation of the pulse and, in addition, there may be some utility in correlating at other lags as well. Consequently, 
a set of transmission code sequences may be stored so that one among the set may be selected based on altitude or 
cell resolution requirements and pursuant to an autocorrelation lag that may be the most appropriate to use for a given 
range of speeds. 

These and other objects and features of the present invention will become more fully apparent from the following 
description and appended claims taken in conjunction with the accompanying drawings. 



EP 0 619 024 B1 



Brief Description of the Drawings 

Figure 1 is a schematic representation of a water-going vessel utilizing means for generating multiple pulse sound 
energy into the water; 

5 Figure 2 is a simplified waveform diagram corresponding to the method of sound transmission shown in Figure 1 ; 

Figure 3 is a schematic representation of a water-going vessel utilizing means for generating monopulse sound 
energy into the water; 

Figure 4 is a simplified waveform diagram corresponding to the method of sound transmission shown in Figure 3; 
Figure 5 is a schematic representation of a typical correlation sonar scenario wherein a water volume, followed 
70 by a seabed, is ensonified by a pulse having a temporal autocorrelation peak at a time lag of 

Figure 6 is a plot of the autocorrelation function of one preferred 300 element coded pulse of the present invention 
wherein the X-axis is defined as code lag in elements and the Y-axis is defined as a unitless autocorrelation co- 
efficient; 

Figure 7 is an example plot of the log amplitude response of an echo at various depths indicating the difference 
is between the true seabed reflection and secondary reflections that are distinguished by a bottom tracking portion 

of the present invention; 

Figure 8a is a scatter diagram of an exemplary current profile produced by a preferred correlation sonar embodi- 
ment of the present invention showing the East/West profile plotted as a function of depth; 
Figure 8b is a scatter diagram of an exemplary current profile produced by a preferred correlation sonar embodi- 
20 ment of the present invention showing the North/South profile plotted as a function of depth; 

Figure 8c is a scatter diagram of an exemplary current profile produced by a preferred correlation sonar embodiment 
of the present invention showing the vertical profile pbtted as a function of depth; 

Figure 9a is a side elevational view of one preferred mechanical assembly for a correlation sonar embodiment of 

the present invention; „_ „.„.._ 

25 Figure 9b is a top plan view of the mechanical assembly for the correlation sonar system shown in Figure 9a; 

Figure 10 is a block diagram of the correlation sonar system shown in Figures 9; 

Figure 1 1 is a flow diagram of the computer software for the correlation sonar system as executed in the electronics 
assembly shown in Figure 10; and 

Figure 1 2 is a flow diagram of the computer software for the maximum likelihood function shown in Figure 11 . 

30 

Detailed Description of the Preferred Embodiments 

Reference is now made to the drawings wherein like numerals refer to like parts throughout. The following Detailed 
Description is organized into sections I - V. Section I discusses alternative pulse generation techniques that may be 

35 used in the present invention. Section II presents the presently preferred correlation sonar hardware. Section III reviews 
maximum likelihood estimation principles. Section IV discusses several methods to calculate the model correlation 
matrix as a function of unknown parameters. Section V describes the operation of the presently preferred correlation 
sonar. The following discussion particularly describes one presently preferred correlation sonar; however, it will be 
understood that the present invention includes other types of correlation systems as well.. 

40 ■ 
I. PULSE GENERATION DISCUSSION 

Figure 1 is an idealized correlation sonar scenario showing a water-going vessel 100 afloat in a body of water 102 
having a bottom 104. The vessel 100 employs a multiple pulse, correlation sonar system (not shown), having two 

45 receivers, to track the bottom 104 and determine vessel velocity. The system transmits a first narrow pulse 106, to 
identify a small range spread of scatterers on the bottom 104, and then repeats to generate a second narrow pulse 
108 to achieve the desired correlation in the transmitted signal at a predetermined processing delay. The two pulses 
106, 108 maintain a distance relationship of t c , the pulse repetition interval, times the speed of sound in water C. 
For instance, Figure 2 demonstrates the concept of a repeated pulse transmission wherein, at two receiver locations 

so on the forward ("fwd") and aft portions of the vessel 100 (Figure 1) : a pulse echo pair 106' and 108", corresponding 
respectively to the first and second pulses 106, 108, will be highly correlated when the ship is moving forward with 
velocity approximately equal to the distance of separation between the two receivers divided by twice the time interval 
between the two transmitted pulses. Note on the top line of Figure 2, time increases to the left (i.e., pulse 108 is 
generated after pulse 106); whereas in lines two and three, time increases to the right (i.e., reception of 106' and 108" 

55 occur after reception of 106' and 106", respectively). 

The resolution of the differential range to the bottom scatterers 104 which can be simultaneously illuminated is 
limited by the width of the pulses 106,103 (Figure 1). This causes the initial portion of the return echo 106' (and 106") 
to come from a center point on the bottom 104 below the sonar system first and then later portions of return echo 106" 
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(and 106*) comes from scatterers on the bottom 104 that are offset from the point and thus at a greater range from the 
system than the center point After a period of time, the system views the scatterers that are beyond the transmitter 
beamwidth and the signal level fades away as indicated by the ramped shape of the final portion of the received echo 
106* (and 106") in Figure 2. At a later time, a second pulse 108 hits the bottom 104 and the scenario repeats. 

5 Because the received signal level attenuates over time, or in other words, the leading portion of the echo is sub- 

stantially different from the trailing portion, the correlation delay of the multiple pulse system (i.e., the time between 
pulses, t c ) should be equal to or longer than the bottom illumination time (i.e., the pulse width). If it were otherwise, 
the system would estimate velocity at significantly reduced precision caused by interpulse interference. 

The timing of the preferred single pulse, or monopulse, transmission of the present invention is quite different from 

io the multiple pulse system. Figure 3 shows a scenario similar to Figure 1 ; however, in Figure 3 the vessel 100 employs 
a monopulse correlation system. A simplified pulse waveform of the presently preferred system is shown in Figure 4. 
The width of a single transmitted pulse 110 is typically long compared to the range spread to the scatterers on the 
bottom 104 (Figure 3). The single return echo 110\110" received by the forward and aft receivers corresponds roughly 
in length to the transmitted pulse 110 and, hence, the scatterers on the bottom 104, inside the beamwidth, are typically 

is all viewed simultaneously. 

In general, the single pulse is coded, e.g., phase coding, to have a pre-specified autocorrelation peak. The mod- 
ulation, or coding, on the return echo 110" of the aft receiver is shitted in time, or delayed, compared to the echo 110* 
received by the forward ("fwd") receiver. This time shift is indicative of the period of time that the vessel 100 traverses 
a distance equal to one-half the distance between the receivers. Therefore, the correlation delay of the presently pte- 

20 ferred system must always be less than the width of the single pulse which is to be contrasted with the multiple pulse 
technique that is typically greater than or equal to the pulse width. 

Figure 5 shows a typical correlation sonar scenario. The water volume in a beam 112 and the seabed 104, which 
has an "effective* tilt as measured by angles a and 0, is ensonified by a coded pulse created by a relatively wide beam 
acoustic transducer (not shown) of width 9 and center frequency f 0 . The nature of a single, coded pulse is such that 

25 the autocorrelation function of the pulse possesses a temporal autocorrelation peak at a predetermined time lag of x c . 
The returning pressure waves are monitored in time at differing spatial locations. Note that at a given instant in time, 
only scattering centers which reside within a series of ring-like structures, e.g., 114 on the seabed 104 are received. 
As time increases, the radius of these rings will increase. The signal return from one receiver at a time, t, will be highly 
correlated with the return from another receiver at a time t+x c provided that the scatters have moved exactly one-half 

30 the spatial vector difference between the two receivers in a time x c . 

This effect is referred to as "waveform invariance". Equivalent^, "waveform invariance" implies that the time domain 
return from a series of stationary scatterers depends not on the exact location of the transmitter and receiver, but rather 
only on the midpoint of the line segment which joins the transmitter and receiver. The principle of waveform invariant 
geometry assumes only that the vertical distance to the scatterers is large compared to the transducer element sep- 

35 aration. For the correlation sonar, the transducer separations are typically of the same order of magnitude as the 
wavelength of sound, and therefore this approximation applies in almost every conceivable case. Another interesting 
feature of waveform invariance is the complete lack of any dependence of the horizontal velocity on the speed of sound 
in the water. 

The signals received by two different transducers can be cross-correlated with a time lag, t c , corresponding to the 
^0 autocorrelation peak of the transmitted pulse. By utilizing a variety of receiver transducers which have different mutual 
spatial separation vectors, it is possible to search for the vector displacement corresponding to the maximum of the 
correlation function. 

The real situation for the correlation sonar is more complicated. Roughly, the location ol the maximum of the 
amplitude of the correlation function contains information about the "horizontal" velocity, while the corresponding phase 

45 of the correlation function contains information about the "vertical* velocity. However, the maximum likelihood signal 
processing algorithm used by the present invention estimates the velocity through a means which does not seek to 
locate the peak of the magnitude of the correlation function nor determine phase under the peak, but rather incorporates 
all available correlations in their full complex form to determine the most likely velocity vector. Additional complications 
must be addressed concerning signal-to-noise ratio (SNR) values and possible sensor pitch and roll changes between 

so transmit and receive. For example, the angles a and p are related to an "effective" seabed tilt and/or transducer tilt 
due to vessel movement. The angles a and p also contain information regarding the beam patterns of the transmitter 
and receivers as well as the angular scattering strength of the seabed. These additional complications tend to slightly 
modify the a and p angles to "effective" a and p angles which are slightly different from the physical angles which 
geometrically relate the sensor coordinate system to the seabed coordinate system. Hence, one key to a successful 

55 correlation sonar velocity log function of the present invention is in the optimal estimation of the unknown three-axis 
velocity vector from the different complex cross-correlation functions formed from the data collected by the numerous 
receivers. 

Referring now to Figure 6, the use of a coded pulse in the presently preferred correlation system provides auto- 
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correlation properties that allows code subsections to be similar but not necessarily identical Moreover, such encoding 
is not used to achieve improved spatial resolution (it is used to improve noise averaging) and thus the waveform is not 
encoded so as to achieve the same function as a repetitive pulse. In fact, non-identical code sect.ons are sometimes 
deliberately constructed to provide zeroes on either side of the desired autocorrelation peak lag which reduce the 

s corruption of velocity estimates which may be caused by autocorrelation side-lobe biases. 

Pulse coding techniques are capable of placing more energy into the water to be available for scattering while 
maintaining the necessary correlation properties at a time lag v It is possible to take advantage of- this technique to 
increase the energy content of the pulse packet while maintaining the high bandwidth of the system. Placing more 
energy per ping into the water is essential to increase the range performance of the present correlation sonar system. 

w Figure 6 shows one preferred encoded pulse. The primary pulse coding technique presently used involves pseudo- 

* random codes which are now generated by a code finding program. The length of the code is specified along with the 
desired time delay for correlation. The computer program, run off-line from the sonar system, then searches for the 
optimal coded pulse sequence. The program performs many optimization searches and keeps track of the best result. 
The corresponding solution is then written to disk. The presently preferred modulation method is a two phase encoding 

is scheme The carrier phase is modulated to be either 6 or 1 80 degrees for each element. Figure 6 shows an example 
300 element pulse code sequence which results in an autocorrelation side-lobe peak at a 1 01 element delay. Preferably, 
each pulse element represents more than two and less than ten carrier cycles. 

The presently preferred code finding process is described as follows. First, a random code is generated using well- 
known methods The autocorrelation function for the initial random code is calculated. A figure of merit based on the 

20 autocorrelation function is then calculated. The figure of merit is related to the sum of squares of all heights divided by 
the square of the height at desired peak or time delay for correlation. The best values for merit are low, and generally 

near 1. . 

Multiple trials are conducted, wherein each trial begins with a new random code. Within each trial the code may 
be mutated and a new figure of merit, based on the mutated code, is compared with the best figure of merit for that 

25 trial If the code is found to be better than the best for that trial, the code elements are mutated at a slow rate, usually 
one at a time, and a more careful search is conducted around the vicinity of the code. If no improvement is found, the 
number of code elements mutated is slowly increased until a predetermined number of code elements is reached. 
Anytime a new trial champion is found, the process slows down to mutating one element to search the neighborhood 
better. After a prespecified number of mutations, the champion from the current random code is saved, and the pro- 

30 cedure begins all over again by generating a new random code. 

Once a best code is found in each randomly seeded trial, it is possible to attempt to tack on some code elements 
at the end which will force lags on either side of the desired peak to be zero. Sometimes this procedure has no solution 
and the code is thrown away. If it is possible, the code elements are appended to the best code and the figure of merit 
for the newly formed code is compared with the existing overall champion from all trials and, if it is superior, it is saved 

35 as the new overall reigning champion. 

Figure 7 illustrates an exemplary echo signal from a sonar system in water having the seabed 104 (Figure 5) at a 
depth of 34 meters. Subsequent to the transmission of a coded pulse, each receiver of the correlation sonar system 
records echo data. Before any correlations or velocity estimation can occur, the data must be searched to determine 
the distance to the seabed, and hence, the portion of the data to process. The seabed is detected by investigating the 

40 amplitude of the response as a function of time. In general, the return from the water volume will decrease from spherical 
spreading and absorption losses. At some point, the pulse will be reflected from the seabed, resulting in a large increase 
in backscattered strength. 

A matched filter approach is used to determine the bottom range from the reflected signal. The presently preferred 
matched filter method is a modified version of that described in U.S. Patent No. 5,122,990 to Deines, et al., which is 

45 hereby incorporated by reference. 

The example field data of Figure 7 plots the logarithm of the amplitude of the received signal as a function of 
apparent depth derived by multiplying #the speed of sound by time. The water return quickly diminishes at the beginning 
of the data stream due in this case to the use of an operational frequency of 76.8 kHz and the cleanliness of the water. 
At a range of approximately 34 meters, the seabed reflection starts. There is a second "bump" at twice the apparent 

so distance which corresponds to a secondary reflection of the initial transmitted pulse down to the seabed, back up to 
the water-air interface, down to the seabed again where it is reflected back to the receiver. Depending on the transmitted 
power, absorption, and depth of the seabed, tertiary and higherorder multiple reflections can occasionally be observed. 
The width of the return is a function of the length of the transmitted pulse and the range to the seabed. Due to thermal 
noise, ambient noise in the water, and flow noise, the receivers always produce a finite signal corresponding to a noise 

55 floor. 
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II. CORRELATION SONAR HARDWARE DESCRIPTION 



As previously noted, the present invention may be embodied as a correlation sonar system having current profiling 
and bottom tracking functions. Figures 8a, 8b and 8c show velocity measurements in a water profile according to a 

5 preferred current profiling function ol the present invention. Figure 9a shows a side elevational view of the mechanical 
assembly, generally indicated at 202, that houses and protects the transducer electro-mechanicals used to implement 
the presently preferred correlation sonar of the invention. An I/O connector (not shown) is located at the back end of 
the transducer assembly 202 so that data can be transferred to the two other subsystems of the correlation sonar 
system 200 (Figure 10) via a transmission cable 220. 

10 A top plan view of the front end of the transducer assembly 202 is shown in Figure 9b. The embodiment of the 

assembly 202 shown has an operational frequency of 22 kHz. However, the transducer assembly 202 is designed to 
be a modular unit so that other frequencies may be specified using different assembly housings according to the various 
sized transducers used as dictated by the selected frequency. Seven transducer elements are packed together in a 
hexagonal configuration (one transducer in the center of the hexagon). All seven of the hex packed elements (five 

>5 labeled 208 and two labeled 210) are used as transmitter transducers. Two of the hex packed elements 210 and six 
other separate elements, arranged in the assembly according to an optimization of vector differences between trans- 
ducers, are used as receiver transducers 210. Each transducer element 208, 210 is approximately 2" in diameter. 
Other numbers, sizes, spatial separations, and transmitter/receiver assignments of transducers could, of course, be 
configured. 

20 The illustrated transducer assembly 202 is approximately 18" in diameter and 11 " long having a weight of under 

150 lbs. and a current profiling range, (with bin resolution somewhat less than a Doppler system), of around 1,200 
meters, about twice the depth of a conventional 75 kHz Doppler system. By comparison, the bottom tracking range of 
one 75 kHz Doppler system weighing about 350 lbs. and having a length of about 1 meter, is around 1,500 meters. 
However, the bottom tracking range of the presently preferred 22 kHz system is 5,000 meters or 16,000 feet. 

2S A functional block diagram of the preferred correlation sonar system 200 includes the transducer assembly 202 

shown in Figure 9 and an electronics assembly 204 coupled to a host computer 214. Other designs may incorporate 
a microprocessor based system inside the electronics assembly 204 to eliminate the need for the separate host com- 
puter 214. In general, the correlation sonar system is conveyed on a vessel, such as the vessel 100 shown in Figure 
5; however, underwater vehicles or fixed mountings, such as platform or buoy mountings, are also possible. 

30 The transducer assembly 202 comprises the transducing elements 208,210 which convert electrical energy into 

acoustic energy for transmitting code sequences, and convert received acoustic energy into electrical energy for anal- 
ysis. The assembly 202 further comprises a set of tuners 21 6 to tune the transducing elements to operate at the proper 
operational frequency and impedance, and, for the receiving transducer elements, a set of receivers 218 to amplify 
the received signal to the appropriate level for transmission to the electronics assembly 204 through an underwater 

35 cable 220 : currently a 19 mm multiply shielded cable, connected to the I/O connector (not shown) of the transducer 
assembly 202. 

The electronics assembly 204 includes a set of four demodulator boards 222 : wherein each board 222 has sine 
and cosine channels for two of the receiving transducers 210. The demodulator boards 222 convert the received signal 
into baseband (low frequency) digital sequences, thus incorporating analog-to-digital (A/D) circuitry. Presently, the 

to complex data is hardlimited and each sample is one of four binary states contained in two bits, one bit per quadrature 
component of the waveform. The time interval between samples depends on the sampling rate. The sampling rate is 
user selectable and is normally chosen to be some factor times the bandwidth of the transducer elements considered. 
Typically, the sampling rate is 2 to 8 times the bandwidth. For example, in the 22 kHz system, the bandwidth is about 
18% of the frequency and each coded element is six carrier cycles but, due to clipping, a sample is made every 1 .5 

^5 carrier cycles. Bandwidth may be increased to obtain more information but at the expense of reduced range due to a 
corresponding increase in the noise. 

A digital signal processor (DSP) 230, which is preferably a 16-bit Texas Instruments TMS320E15, accumulates 
the sequences. The DSP 230 forms cross-correlations of sequences received from the eight receivers into two or more 
3x8 matrices, i.e., one correlation coefficient per receiver transducer pair, at each of the time lags, typically 0 and x c . 

50 A timing generator 232 controls a power amplifier 226 to generate the appropriate coded pulse, which is generated 

into the water via the transmitting transducers 208. The electronics assembly 204 is powered by a power supply 236 
which is regulated by a power regulator 234. The preferred correlation sonar system 200 requires less than 100 Watts 
of power. 

The operation of certain other functional blocks, primarily the demodulators 222, the DSP 230 and the timing 
55 generator 232, in the electronics assembly 204 is coordinated by a microcomputer or CPU 228, which is preferably a 
CMOS 68000 microprocessor available from a number of vendors including Motorola. The host computer 214 com- 
municates with the electronics assembly CPU 228 across a line 238 : presently an RS-232 serial communications line, 
via an I/O port 240. The I/O port, or serial I/O board 240 connects to a computer or PC 242. It is also possible to replace 
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the CPU 228 with an interface card for the electronics assembty. This interface card is then connected via a cable to 
a PC interface board located in the host computer 214. Communication between the host computer 214 and various 
electronics modules (e.g., DSP 230. timing generator 232, demodulators 222) is then P^jble without the » use > of a 
CPU card 228 The computer 242 also includes a secondary data storage media 244, preferably eithera 100 Megabyte 
hard disk drive, a 44 Megabyte Bernoulli box or an optical disk. The computer 242 preferably includes a floating point 
processor and is presently configured with a Number Smasher-860 board, available from Microway of Kingston, Mas- 
sachusetts, which includes an Intel i860 64-bit, RISC architecture, microprocessor. 

The computer 214 provides a convenient interface between the user and the correlate sonar system 200 by 
controlling various aspects of transmit, receK/e, storage, and analysis of data. The computer 214 allows the user to 
select a transmit pulse from a series of coded sequences stored on the data storage device 244. as well as define the 
.ength of the receive signal to sample. The receK/ed signal is then divided into range gated bins where.n each bin 
contains a water velocity which has been estimated by maximum likelihood signal processing which is described below 

III. MAXIMUM LIKELIHOOD SIGNAL PROCESSING 

Subsequent to the transmission of a phase modulated pulse, each of the receivers 210 (Figure 10) record inbound 
pressure waves (complex values of voltage) as a function of time. The data is divided into bins wherem each bin is 
used to calculate a velocity estimate. The method outlined in this section is considered to be excellent at estimating 
the velocity: however, it is very computationally demanding. 
20 m all the analysis presented here it is assumed that the incoming pressure is a complex Gaussian random variable 
as a function of time This is justified because of the coupling between the central limit theorem and the large number 
of randomly placed scattering centers in the range bin of which the water velocity is being estimated. The general 
expression for the probability density function (PDF) of observing the data vector z given the parameters were in the 
state, 9 for n complex data points is 



2S 

p(z;6) = J- n jl- | exp(-zVz), 0) 

71 

30 where R is the n X n correlation matrix describing how data point 3 is correlated with data point ^ IRI is the determinant 
of R and p{z- 6) is the probability of obtaining 2 given that the parameters were in state 0. With a prion assumptions 
concerning the probability distribution for the parameters, the most likely set of parameters can be obtained through 
solving for 6 in 



WlD-o^MIl^) (2) 



where p(9:z) is the probability that the parameters were in state e given the data was observed to be z. 

40 A closed form analytical expression for the form of the cross-correlation function for water column returns enables 

the model correlation matrix, R, to be compuied as a function of the unknown parameters (including velocity) . This R 
matrix may then be used along with the data correlation matrix, Zj 2^, to form the probability function. The values of 
the unknown parameters may then be varied in an iterative search for the maximum of the probability function, which 
is equivalent to solving the set of simultaneous equations represented by equation (2). The expansion involves infinite 

45 sums over hall-integer spherical Bessel functions and Legendre polynomials. In practice the sums may be truncated 
after several terms to obtain an adequate approximation. 

To simplify the calculation, the natural logarithm is applied to equation (1 ) and. as a minimum is desired, the equation 
is multiplied by -1 so that the function of unknown parameters, g, becomes, 



9( v x> v r v r SNR i=Q> SNR T=x c ' a $' W,a)=Anp(z;Q)= 

_l n JL + ln| + 77ace( R^Ff^ ( 3 > 

The first term can be ignored in the calculation as it is an offset curve. In one presently preferred system 200, the set 
of unknown parameters are as follows: 
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v x = the water velocity along the x-axis; 

v y = the water velocity along the y-axts; 

v 2 = the water velocity along the z-axis; 

SNR^ = the signal-to-noise ratio at lag x=0; 

SNR T=XC = the signal-to-noise ratio at lag x=x c ; 

a = effective relative angle of the azimuth (Figure 5); 

P = effective relative seabed tilt from the vertical; 

W = width parameter of angular scattering function; and 

a = wing parameter of angular scattering function. 

Thus, in the maximum likelihood method, the unknown parameters are varied to minimize the function, g, and thereby 
maximize the probability of observing the received data in a specified parameter state. 

IV. MODEL CORRELATION MATRIX 

The model correlation matrix Rmodei- or J ust R in tnis section, can be calculated using a closed form expression 
for the spacial correlation function R. As water volume returns and bottom returns require somewhat different processing 
due to seabed tilt and seabed roughness issues, the following discussion presents an expression for the spacial cor- 
relation function used with water volume returns followed by several expressions for the spacial correlation function 
used with bottom returns. The expressions used with bottom returns differ as folbws: the first expression assumes the 
seabed exhibits a Lambertian scattering function, i.e., cos(0); the second expression assumes the seabed exhibits a 
Kirchhoff approximation scattering function, i.e., a(9); and the third expression is a heuristic which much simplifies the 
calculations for R. 

In the presently preferred embodiment, the structure of R is determined in the following way. The inbound pressure 
data points in time are numbered 1-8 for data at time t=0 for receivers 1 -8, respectively, and 9-16 for data at time t=x c 
for receivers 1 -8, respectively; then 17-24 for data at time t=At (where At is the sample interval) and 25-32 for data at 
time t=r c + At, again for receivers 1-8, respectively, and so forth. A very large correlation matrix of size 8m x 8m is then 
generated for a system comprising 8 receivers and collecting B m' data points in time at a sampling interval of At. The 
value for "m" depends on a variety of factors which may include altitude, water profiling bin resolution, and number of 
code elements transmitted. For simplicity, it will be assumed that the sampling interval At is commensurate with the 
system bandwidth so that each sample in time is independent. It is then easily seen that the 8m x 8m large matrix is 
really a block diagonal matrix with 16x16 submatricies down the main diagonal and zeros elsewhere. That is, data 
points 1-8 are correlated among themselves (spatial correlation function at 0 timeshift) and they are correlated with 
entries 9-16 (spatial correlation function at x c timeshift) but they are uncorrected with all other data points. Likewise, 
entries 17-32 constitute a 16 x 16 submatrix, and so forth. Because the monopulse approach typically ensonifies all 
the scatterers in the beam pattern simultaneously, the nature of each 16 x 16 submatrix along the main diagonal is the 
same, i.e., each 16 x 16 submatrix along the main diagonal of the large 8m x 8m matrix represents the spatial correlation 
function evolving in time. However, an advantage of the beam filled monopulse method is that the spatial correlation 
function becomes stationary, i.e., not a function of time. Thus, only an individual 16x16 submatrix need be considered 
for R. The spatial correlation function at zero time shift is contained in the upper left 8 x 8 and lower right 8x8 portion 
of R. The spatial correlation function at x c time shift is contained in the upper right 8x8 and lower left 8x8 portion of 
R. R is of course, Hemnitian. 

As is well-known in the transducer technology, the beam pattern, W(0), where 8 is the angle of measured energy 
away from the axis of wave propogation, can be represented as a finite series of Legendre polynomials with coefficients 
b n , which thus allows the spacial correlation function to be a function of sum of integrals which involve the n th term of 
the beam pattern. (One description of Legendre polynomials is provided in Mathematical Methods for Physicists by 
George Arfen, Academic Press, New York, pp. 534-608.) Calling each of these integrals R n , 




J, 



P k (cos6) 



X/2*k 



(6) P k (cose ) 



(4) 



where J is a Bessel function and P is a Legendre polynomial. Thus, the spacial correlation function is the sum of 
integrals as follows: 



(5) 
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R = £ A, *„ 

5 

where each b n is a Legendre polynomial coefficient, and where 

y = k(v x (a4-a 1 )+v y (a 5 -a 2 )+v z (a 6 -a 3 )J/c ) lhe expression for a correction to the phase corresponding to path differ- 
ences of order v/c, where v is water velocity and c is the speed of sound in water 
w k = 2vJX 

c = speed of sound in water 

X = wavelength of the propogated sound wave 

and further where the first and second receivers are located in the sensorframe at (A, ,A 2 ,0) and (A4, As.O) respectively. 
15 Because the sensor can change attitude between the time of transmission and the time of reception, a rotation matrix 
must be applied to obtain the receiver locations appropriate in the coordinate system of the sensor at the time of 
transmission: 

2o ' a^ =cosct*cosp. [A 4 -A l )cos7-(A 5 -A 2 )siny]-sina[(A 4 -A 1 )sinY+(A 5 -A 2 )cosY] 

a s -a 2 = sina-cosPtCA^j-A^cosy-tAs-AgJsiny] + cosa[(A 4 -A 1 Jsiny+fA^A^osy] 

a 6 -a 3 = -sinp[(A 4 -A 1 )cosy-(A 5 -A 2 )siny] 

25 ^ 

a -tan (-siny-cosp/cosy) 

222 1/2 
5 = Al(2TV z +a 6 -a 3 ) + (2*^+3^) + (2w y +a 5 -a 2 ) ] (6) 

30 2 

[(2t v-a 4 -a,) +(2t V a 5 -a 2 ) ] 

tane = ^ \') 

2xv x +a G *a 3 

Unfortunately, the Legendre polynomials are orthogonal on the interval (-n/2,7c/2), and not the integration interval of 
35 (0,7t/2) . However, after some further mathematical work on the equation it can be shown that the closed form expression 
for FL is as follows: 



R - = e 



40 



45 



6 *" mR% ■ 222*1 

[ll]^ 2 £ (1/2**) z- k J., 2 . k l6) 

00a 

x P^cose) *( ! . _ if n is even 

n! ! (>c-l) ! ! Ik (k+1 ) -n (n + 1) J 

[•2JTJ1/* £ (i/2*Jc) i-* c7 1/2 . k (6) 

evtn v 

x P k [cose) (-1) <»•*- l,/2 "_— , JiVtM 1 ?! ! wu ,m if n is odd ) 



where 

5 = At(2t^ + a 6 -a 3 ) z + (2tv x+ a 4 -a,) z + (2xv y+ a s -a 2 ) z J (9) 
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cose - 



2x^+a 6 -a 3 



(5"0 



(10) 



10 



Depending on the parity of the beam pattern term being considered, i.e., n is even or odd, only one of the sums 
in equation (8) is to be used. 6 is the total distance away from the bistatic point in all three dimensions and e is a phase 
angle describing the relative contribution of vertical and horizontal displacement. In practice, the sums in equation (8) 
may be truncated after a few terms to obtain an adequate approximation to F^. Equation (5) is then used to sum over 
the coefficients of the beam pattern to obtain the final, total correlation value. 

A spacial correlation function to model bottom returns coming from a seabed is also important however such a 
function should incorporate the existence of a tilted seabed as, for example, shown in Figure 5. A four variable closed 
form expression is provided below: 



is 



20 



25 



where 



R « 



X [ 



2£ b ± i; o P 1 (cos0) 

2-0 

(27i/r) l ' 2 £ Jbj [P^cosfi) £ i-J I,\ J JmU2 lr) P^cosij)* 

2£ Pf(cos/3) (-1)" cosfm(0 o -cr)] £ i ~ } J" * «/j. l/2 (•*") P/(costj)]] 

(11) 



Variable 



Significance 
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r 

tanq 



P 



-a 



(2ic/A,) X the total three axis distance the lag is from the bistatic point (2tv x: 2xv y , 2tv z ) 

the ratio of the horizontal distance away from the bistatic point to the negative of the vertical lag away 

from the bistatic point (defined in the seabed frame) 

the nadir tilt angle of the seabed 

the azimuthal angle corresponding to the ratio of distance away in the y lag to x lag indexed off of the 
azimuthal angle in which the seabed is tilted 



and where ring geometry coefficients are defined as: 



If, = <l/2+j> \\ m \ ; ) j ' m} ; f sxnd cosO Pf (cosS) Pf (cos*) dd 

(12) 

In equation (12), where 1 and j are dummy indices which are being summed over, it is assumed that the seabed 
scattering follows a Lambertian law (i.e., intensity is proportional to cos8) and that the function includes the combined 
effect of a series of rings from 9 1 to 9 2 . If the angular scattering law of the seabed is not cos e, but a(8) instead, then 
a{6) must be substituted for cosG in equation (12). The expression of equation (11 ) may be too complicated for many 
real-time applications so a two variable approximation has been derived as follows: 



(2ir/r) 1/2 t*** l m 

■2Z b.I° 10 J '° v0 (12.1) 

ss J ' e 



Unfortunately, under certain seabed roughness conditions the scattering function may not follow a Lambertian law. 
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This may be true it the scattering function of the seabed is strongly dependent on the angle of the beam, being strongest 
at small angles near nadir Thus, some alternative to the Lambertian function may be required to more accurately 
predict acoustic scattering on a seabed. 

The Kirchhoff approximation assumes that the seabed surface is slowly varying so that the radius of curvature is 

5 much greater than a wavelength. Some researchers, including Clay, C.S. and Medwin, H. 1977, Acoustical Oceanog- 
raphy, John Wiley and Sons Publishing, have noted that the Kirchhoff approximation works reasonably well under a 
variety of realistic situations provided the incident angles do not approach grazing angles. Fortunately, for almost any 
conceivable correlation sonar geometry, most of the sound is closer to normal incidence than grazing incidence. 
The expression for the backscatter cross section per unit area per unit solid angle for the diffuse scatter field under 

io the Kirchhoff approximation can be derived, (following the development by Jackson, et al. in "Application of the Com- 
posite Roughness Model to High-Frequency Bottom Backscattering", J. Acoust. Soc. Am. 79 (5), May 1986) and is 
defined as follows: 



o(6) - £0 fw^fwy eX pf2iJc sin(fl> x] (13) 

27T 2 COS 2 0 J o I 

x exp[-2Jc 2 cos 2 (fl) D(r)] 

20 where the structure function D(r) is used instead of the correlation function and the structure function is defined as D 
(r) = 2[B(0) - B(r)] where B(r) is the spatial autocorrelation function of the seabed and is a reflection coefficient. 

Most data on bottom roughness are derived from bathymetry and have a resolution of 1 00 meters or larger. Two- 
dimensional roughness spectra with centimeter scale resolution have be measured by Akal, T and Hovem, J., 1978, 
"Two-dimensional space series analysis for seafloor roughness," Marine Geotechnol. 3, 171-182. Others in the litera- 
ls ture have provided one-dimensional roughness spectra with centimeter-scale resolution. The data from course and 
fine seabed roughness measurements agree in having power-law roughness spectra. Therefore, this model assumes 
that the two-dimensional seabed roughness statistics are Gaussian and isotropic with a spatial power spectrum of the 
form 



<D(*)=P* 0 4 ) 



where k is a two-dimensional wave vector having a magnitude equal to the wavenumber k, and p and y parameters 
describe a particular power spectrum of the seabed. The p and y used in equation (14) are different and should not be 
35 confused with the p and y angles used earlier. The noted high resolution data has the following approximate range for 

y:3<y<3.5. 

After relating the structure function to the spatial power spectrum and the substitution of a = (y/2)-1 is made then 
4 o ' D(r)=cf h r 2ci (15) 

where 

45 ^> = 2?tp 2 2a T(2-a) 

h a(1-cc)r<1-a) 

a = (y/2)-1 (15.1) 

50 where C h is a structure constant, also sometimes called a structure function parameter, and T(x) is the gamma function. 
Again a here should not be confused with the a angle used earlier. The final backscatter equation can be derived by 
incorporating the structure function given in equation (15) into equation (13) and simplifying to yield 

55 
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R 2 

o(d) = If f du u J c (u) exD(-gu 2a ) (16) 

47T sin 2 (d) cas 2 {6) J 0 ^ 



where 

g = 2°" 2a) ** 2 " 2a) cos 2 (6) sin- 2a (6) (17) 



where q is an arbitrary variable chosen to simplify the equation. Equations (16) and (17) can be numerically computed 
to yield the backscatter cross section per unit area per solid angle as a function of angle 0. 

Although the Kirchhoff approximation appears to provide an accurate estimate of scattering, in investigating the 
J 5 various scattering curves corresponding to different degrees of seabed roughness a heuristic for the spacial correlation 
function R was discovered. The heuristic is used as an alternative expression when there are certain real-time con- 
straints in the processing. The heuristic is a general function for band shapes as can be found, for example, in Spectral 
Analysis by James Blackburn. The magnitude of the complex R, the spacial correlation function, expression is as 
follows: 

20 

R= j 2 < 18 ) 

(1 +(2* -vi2(x-x Q yw\Y a 

25 

where Xq is the center of the curve, W is the width parameter defined as twice the distance from the center to the 1/2 
curve height point and a is a wing parameter specifying the extent of the curve wing. R peak is either initially set high 
to .9 or estimated from the amplitude pressure data of the seabed return relative to the noise floor and a is set to 1/2. 
Note that if a is set to 0, the curve is a Gaussian, and if a is set to 1, the curve is a Lorentzian. In addition to the 
so magnitude of the complex R shown in equation (18), the heuristic includes a phase angle factor which is applied to the 
magnitude to generate the complete heuristic complex R. The phase angle factor is as follows: 



35 



V. OPERATION OF THE CORRELATION SONAR 



Maximum likelihood processing specific to the presently preferred acoustic correlation current profiler (correlation 
sonar) embodiment of the correlation sonar system is presented by general reference to Figures 11 and 12. In the 
Figures, function names used in the preferred software appear to the right of their corresponding state. Software for 
the present system is written in the FORTRAN language. The software described herein, particularly the portion listed 
in the attached Microfiche Appendix, was translated from source code to machine-readable object code using the NDP 
FORTRAN Version 4.0b compiler licensed from Microway. Nonetheless, one skilled in the technology will recognize 
that the steps in the accompanying flow diagrams can be implemented by using any of a number of different computer 
languages and language translators. 

Now referring more specifically to Figure 11, there is shown the correlation sonar software, generally indicated at 
250, that is stored in the memory (not shown) of the host computer 21 4 (Figure 1 0), and which controls signal generation 
and echo sampling, and conducts analysis of the sampled data to obtain velocity measurements. Beginning at power 
on, the correlation sonar enters a start state 252 and continues to state 254 to initialize the software by creating data 
structures and reading parameters from files on the disk 244, and at state 256 data is stored in registers to initialize 
the hardware. 

At present, one of the software initialization parameters specified by the user is "bottom track mode". If the param- 
eter is true a bottom track has been requested, and the execution flow proceeds down the right portion of the flow 
diagram; otherwise, a water profile has been requested, and the execution flow proceeds down the left portion of the 
flow diagram. Of course, one skilled in the technology will comprehend that the modes are for programming convenience 
and that the functionality of both modes could be combined so that bottom tracking and current profiling occur 'together". 
In that way, the earth reference velocity of the water currents could be completely obtained within the correlation sonar 
200 by subtracting the bottom tracking, or vessel, velocity from the current velocities. 
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As the right and left portions of the flow diagram in Figure 11 overlap in most functions, only the left side will be 
discussed herein on a state-by-state basis. Thus, assuming that the user has specified measurements for a water 
profile, such as the example shown in Figure 8, the correlation sonar enters state 260 to select a specific code that is 
used to control the timing generator 232 (Figure 232) and, hence, generate a coded pulse from the transmitter trans- 
ducer 208. 

The sequence data indicative of the echo signal is received from the demodulators 222 (Figure 10) by the CPU 
228, and forwarded to the host computer 21 4. The correlation sonar 200 then enters state 264 to unpack the sequence 
data. Specifically : the sine and cosine binary data for two receiver transducers 210 and 1 byte of amplitude data, 
sampled at a lesser rate, for each transducer 21 0 are all packed into the digital sequence. Therefore, the different data 
and data from different sources is extracted at state 264. Proceeding to state 266, the cross^correlation matrix, R(0) 
for a time lag of 0 is created for all eight of the receiving transducers 210 (Figures 9 and 10), thus creating an 8 x 8 
matrix, wherein each row and column corresponds to a receiving transducer. Subsequently, at state 268, the cross- 
correlation matrix, R(x) for a time lag o1 i c is created as another 8x8 matrix. Then, a "seed" is generated for the 
maximum likelihood signal processing at state 270. Here, accuracy is not the consideration so much as providing an 
initial guess for the water velocity unknowns, v x , v y , and v 2 . This is accomplished by using a 3 x 3 matched filter approach 
on the R(t) matrix to average the contents of every nine neighborhood amplitude correlation coefficients. A first moment 
technique is then used and the results are divided by 2r to obtain an initial estimate of v x and v y . The velocity v z could 
be set to zero. However, more preferably, the phase shift of the highest correlation coefficient can be used to estimate v z . 

Two transducer pairs, one up from the origin and one to the right of the origin, are used to index the R(0) matrix 
so as to estimate effective sensor tilt or seabed angles a and p. The phase shifts are indicative of the direction of sound 
propagation in two dimensions. 

The initial SNR 0 estimate can be obtained from the amplitude pressure data relation to the noise floor, or it can be 
set to a high value, say 10, arbitrarily. The SNR^ estimate can be determined by selecting the highest correlation 
coefficient from the R(t) matrix and then using the following. equation;. . . . . 

] _i . 1 n g\ 

magnitude ( p ) SNR v 1 

After setting the above-noted estimated parameters, the correlation sonar 200 enters state 272 to handle maximum 
likelihood (ML) signal processing. The inputs to state 272 can categorized into three groups, group one contains the 
following inputs: a, p, SNRq, SNR t , v x , v y , and v z; (w and a are also included for bottom tracking), which were calculated 
at state 270. Group two contains the changes in the group one parameters. Initially, these values are determined based 
upon reasonable guesses. Group three contains cross-correlation matrices R(0) and R(t). The outputs of the ML state 
272 are the new group one parameters. The ML processing is executed on the host computer 214 (Figure 10) after 
either receiving the raw data and forming the correlation matrices or receiving the cross-correlation matrices directly 
from the DSP 230 of the electronics assembly 204. The function of state 272 will be discussed hereiribelow with ref- 
erence to Figure 12. 

Continuing to refer to Figure 11 , after the completion of the ML signal processing (state 272) the correlation sonar 
200 moves to state 274 to store the outputs: a, p, SNRq, SNR^, v x , v y , and v z , (and w and a for bottom tracking), to the 
data storage device 244 (Figure 10). Also, -ln(P) is stored which is indicative of a goodness of fit. Then, at state 276 
a coded pulse is selected, which is not necessarily the default code as was selected at state 260. In fact, the code is 
selected according to the relative measured velocity of the water and in the case of bottom tracking, the currently 
estimated altitude to the seabed. 

The pseudorandom codes are selected primarily based on the length and correlation lag time. Length is important 
for the following two tracking criteria: 

(i) for water tracking because it determines the spatial resolution of each water bin. Long codes will produce coarse 
resolution in the profiling of current velocity. They, however, have the advantage of depositing the most amount of 
energy into the water and permit maximum range. Shorter codes allow better resolution (up to a point where the 
beamwidth geometry begins to dominate the bin profile resolution), but they do not place as much energy into the 
water. 

(ii) for bottom tracking the code length is chosen to be long enough to typically ensonify all the scatterers within 
the beam pattern on the seabed and allow sufficient time to perform the necessary cross-correlations. 

Lag time : i c , is important because the location of the peak is approximately given by 2x c velocity. Thus, x c is chosen 
based on velocity such that the peak location will remain on the configured array of receivers, otherwise accuracy in 
the estimated velocity is compromised. However, if t c is selected to be too low, thereby virtually guaranteeing the peak 
is on the array for even high velocity precision in the estimate tends to be lost. These two considerations must be 
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constantly monitored and traded-off to arrive at a proper selection of lag time. 

The ping for the selected code is generated at state 278 in the same manner as was described with respect to 
state 262. In receiving the echo signal for current profiling, the echo return signal is range gated (or, equivalently, time 
gated) by looping over the states 280-294, wherein state 280 is a counter state for incrementing the bin, or cell, number. 
The bin number is calculated according to the length of the pulse and the velocity resolution specified by the user The 
loop is finally exited at state 294 when a test of the bin counter indicates that all bins for this ping have been processed. 
It should be noted that preferably a weighing function based on depth is used to determine the number of profile bins. 
For example, the program DETPRO in the Microfiche Appendix will provide a weighted bin profile. Velocity measure- 
ments based on a new ping are initiated again at state 276, and they continue until the correlation sonar 200 is either 
reset or powered down. 

Returning, in the discussion of Figure 11 , to the right portion of the flow diagram which relates to the bottom tracking 
mode of operation, most states have been labeled with the same number as the left portion (current profiling mode) 
followed by prime f) to indicate identical or similar functionality. However, states not shown in the left portion include 
states 296, 298 and 300. Specifically, state 296 is entered after an initial ping (state 262 1 ) to find the seabed. If the 
seabed is not found, a new code is selected and ping generated until it is found. This state is accomplished using a 
matched filter approach which is a modified version of the process disclosed in the common assignee's previously 
cited patent. The loop from states 276*-292' is continuously executed unless the seabed is not found, e.g., the vessel 
passes over a canyon. In that case, the correlation sonar 200 moves from state 298 to state 300 to report the loss of 
tracking and a the bottom track initialization (states 260',262\296) is again carried out. 

One presently preferred function for maximum likelihood (ML) signal processing is shown diagrammatically in 
Figure 12. Note that parameters characterizing seabed roughness, e.g., the width parameter, W, and the wing param- 
eter, a, both of which were discussed above, are not shown in Figure 12, but could be readily incorporated therein by 
one who is skilled in the relevant technology. 

Generally, the function uses a dual fit approach, namely, three parameters for water profiling or five parameters 
for bottom tracking are fit, then the remaining four parameters are fit. In particular, the effective tilt angles, a and p, 
and the signal-to-noise ratio at time lag O, SNRq naturally reside in the lag 0 correlation space. Likewise, SNR t , and 
the water velocities, v x , v y , and v 2 , naturally reside in the lag x correlation space. 

The ML function, used in states 272,290,272', 290' (Figure 11 ), is entered by the correlation sonar 200 (Figure 1 0) 
at start state 310. Whether the first or second parameter pass is being calculated depends on the FITFLG (fit flag) 
which is initially set to true at state 312 to indicate a, p and SNR 0 will be fit. As the multi-variate, non-linear programming 
optimization relies on the well-known simplex algorithm, a four (six for bottom tracking) vertex polyhedron is formed 
by evaluating the negative of the natural logarithm of the probability function, g, , at the model parameters and at each 
parameter varied by an offset as follows (brackets denote additional parameters required for bottom tracking) : 



g,(a 0 ,$ 0 ,SNR 0 ,[w 0 ,a 0 }) (20.1) 

g A (a Qt p Q + AP 0 , SNR Q , [w 0 , a Q }) (20.3) 

g A (a Q: p 0 , SNR Q + ASNR Q , [w Q , a Q )) (20.4) 

g, (a 0 + Ao, P 0 , SNR Qt [w 0 , a Q ]) (20.2) 

(a 0 , p 0 , SNR 0 , [w 0 + Aw Qt a 0 )) (20.5) 

(a 0 , p 0 , SNR Q , [w 0 , a Q + Aa Q )) (20.6) 



It is important to note that the offsets are estimated. If the offsets are too large, the simplex method, which iteratively 
moves along the edges of the polyhedron, may be misled into a local minimum instead of the desired global minimum. 
On the other hand, if the offsets are too small, the simplex method may take too long to complete because the starting 
point of the iteration is far from the global optimum and the incremental change at each iteration is so small. 

The simplex method in the correlation sonar 200 is presently implemented in a FORTRAN subroutine named 
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AMOEBA (since the method "oozes" down the edges of the tetrahedron). The subroutine is based on the simplex 
method given in "Numerical Recipes, The Art of Scientific Computing" by Press, Flannery, Teukolsky, and Vetteriing, 
specifically described in Section 10.4 entitled "Downhill Simplex Method in Multidimensions"; Press Syndicate of the 
University of Cambridge 1 . One of the parameters in the subroutine is ITER, the number of iterations to be performed 
5 by AMOEBA. At each iteration the unknown parameters a, p, and SNRq (and [w c , aj for bottom tracking) are evaluated 
in g-, as follows: 

£>, =1" Id* +Tia«(* w £ >dW ) (21) 

10 

where A is an 8 x 8 matrix equal in structure to R(0). Thus, a user defined maximum number of iterations will prevent 
an endless iteration over the search space, and a test is made at state 31 8, to determine whether the maximum number 
of iterations was exceeded. If so, a report of the status is written to the video display of the host computer 214 (Figure 
10) and the ML process completes with incomplete results at state 322. 
is Otherwise, if the simplex method (state 316) did converge within the predetermined number of iterations, the 

presumed global optimum forthe set of unknown parameters in the first gamma function, g, , is stored to system memory 
at state 324. Then, the FITFLG is set to false, at state 326, to indicate that the second set of unknown parameters will 
be fit. In this instance, at state 328 the unknown parameters, SNR,. , v x , v y , and v 2 are fit according to the second 
gamma function, g 2 , as follows: 

20 

ff 2 = ln|cW(P mBd J|+7>ace(P <fcl ,p- ( l 0lW ) (22) 

A five vertex simplex is created by evaluating at the five. parameter sets defined r by._varying the.pffseis. similar to 

25 the method employed in state 314 described above. In calculating g^, a 16 x 16 matrix, P (or RHO), is created wherein 
the matrix can be described as having four 8 x 8 complex submatrices. The upper left submatrix is R(0), the upper 
right submatrix is R(t), the lower left submatrix is the complex conjugate of R(t), and the lower right submatrix is R(0). 
This special matrix structure allows a fast matrix inversion. 

Again, the simplex method is applied to search a multi-variate parameter space, this time a five parameter space 
30 at state 330. If the maximum number of iterations are determined to have been exceeded at state 332, then a report 
is made (state 334) and the correlation sonar 200 exits the ML procedure at state 336. Alternatively, continuing from 
state 332 to state 338 a solution has been found and the control flow continues to return the seven (nine for bottom 
tracking) parameters in system memory. The ML process then returns control flow back to the correlation sonar top- 
level flow as shown diagrammatical ly in Figure 11. 

35 

Claims 

1. A correlation sonar system (200) having a transmitting transducer (208) for emitting acoustic energy onto a re- 
40 fleeting surface, a pulse generation means for generating a single pulse (110) of electrical energy to the transmitting 

transducer a plurality of receiving transducers (210) for receiving a pulse echo (110') from the reflecting surface, 
a plurality ol demodulators (222) : each demodulator connected to one of the receiving transducers and generating 
a set of values and maximum likelihood means for estimating velocity using correlated values, characterised by 
further including: 

45 a means (230) for auto-correlating said values for each receiving transducer, and cross-correlating between 

different receiving transducers, at one or more predetermined time delays less than the pulse width. 

2. The system of Claim 1 , wherein the pulse comprises a plurality of coded elements. 

50 3. The system of Claim 2, wherein each coded element comprises a plurality of carrier cycles. 

4. The system of Claim 3, wherein each coded element is more than two and less than ten carrier cycles. 

5. The system of Claim 2, wherein the coded elements are phase encoded. 

55 

6. The system of Claim 5, wherein the coded elements are encoded with either of 0 degree and 180 degree phases. 
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7. The system of Claim 2, wherein the pulse includes sets of non-identical code elements so as to provide zeroes 
on either side o1 a preselected autocorrelation peak lag. 

8. The system of Claim 1 , wherein the pulse has a temporal autocorrelation peak at the predetermined time delay 

5 

9. A correlation sonar system (200) having a means for generating a pulse (110) into a liquid having reflective surfaces, 
a plurality of transducers (210) for receiving the echo signal (110') created by the pulse contacting the reflective 
surfaces, and maximum likelihood means for maximizing the probability of a set of preselected parameters, in- 
cluding at least one component of velocity, having sets of correlations indicative of a parameter state, said system 

10 characterised by further comprising: 

correlation means for forming correlations of the echo signal received by each transducer at a predetermined 
time lag with the echo signal received by the same transducer at one or more other time lags, and with the echo 
signals received by the other transducers at one or more of the other time lags, said correlations providing for the 
spatial relationship of the transducers. 

75 

10. The system defined in Claim 9, wherein the correlation means includes means for calculating correlations of the 
echo signal using a spacial correlation function defined as follows: 



ji-0 



25 where 

b n is a Legendre polynomial coefficient. 

R n is a spacial correlation function involving the nth term of the beam pattern W(6), 
W(9) is the system beam pattern, and 
30 e is the angle of measured energy away from the axis of wave propagation. 

11. The system defined in Claim 10, wherein the closed form analytical expression for R n is defined as follows: 



R n = e^{[i^] 1/2 (l/2+n)i~ n J l/2mn i6) PJcose) —1 



27T - 



[11]^ £ (1/2+*) i-* 



k*l 

x P,(cose) (-1) <-*-»/2 *l 1 ; r-rr if ™ ^ even 

nl ! (k'l) I ! [k(k+l) -n in+1) ] 

£ (1/2+k) i'* J ua . k l6) 

U Jc-0 

45 even 

x P fc (cose) (-D ""*- 1 >/ 2 . {k ~] ] 8 1 n [ ! if n is odd } 



so where 

n is the nth term of the beam pattern, 

\j/ = k[v x (a 4 -a 1 )+v y (a 5 -a 2 )+v 2 (a 6 -a 3 )]/c,the expression for a correction to the phase corresponding to path dif- 
ferences of order v/c, where v is water velocity and c is the speed of sound in water, 
55 s - (2tUX) x the total three axis distance the lag is from the bistatic point, where A. is the wavelength of the 

propagated sound wave, 
J is a Bessel function, 
p is a Legendre polynomial, 
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e is a phase angle describing the relative contribution of vertical and horizontal displacement, and 
k = 2k/a, where X is the wavelength of the propagated sound wave. 

12. The system defined in Claim 9, wherein the closed form analytical expression for the spacial correlation function 
includes a Lambertian scattering law function. 

13. The system defined in Claim 12, wherein the spatial correlation function is defined as follows: 

R = f x [ 

1*0 

(27r/r) 1/2 £ Jb 2 [^(cos/J) £ Ii, Jj. U2 (r) P,(cosij) + 

i«0 J*Q 

1 • 
2^ Pf(cos£) (-1)- cos [m(<t> 0 -a)) £ i"J J^ 1/a <r) P/tcos??)]] 



where 



\jz = = k[v x (a 4 -a 1 )+v y (a 5 -a2)+v 2 (a 6 -a 3 )]/c,the expression for a correction to the phase corresponding to pat 

differences of order v/c. where v is water velocit and c.is the speed of sound in water, 

L is an upper limit to a summation, 

b n is a Legendre polynomial coefficient, 

l]j ring geometry coefficients defined as follows: 



iTi = d/2+j) \\' m \ ! { A' m \ ! f sin6 cosB Pf(cosfl) Pf (cosfl) d6 , 



where 



I is a dummy index which is being summed over, 

j is a dummy index being summed over 

m is a dummy index which is being summed over, 



P is a Legendre polynomial, 

P is the nadir tilt angle of the seabed, 

r = (2n/X) x the total three axis distance the lag is from the bistatic point (2tv x , 2xv y , 2tv z ), where X is the 

wavelength of the propagated sound wave : and 

t is a time shift, 

J is a Bessel function, 

cosr| is the negative of the vertical lag away from the bistatic point. 

<t> 0 -a is the azimuthal angle corresponding to the ratio of distance away in the y lag to x lag indexed off of the 
azimuthal angle in which the seabed is tilted. 

14. The system defined in Claim 12, wherein the spatial correlation function is defined as follows: 



R~ J-zJLL—. x [I^Ei-^J.^tr) P^cost;)] 

2Eb i J J a . 

1*0 
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where 



r is (2n/X) x the total three axis distance the lag is from the bistatic point (2tv x , 2w yi 2tv 2 ), where X is the 
wavelength of the propagated sound wave and 

t is a time shift, 

v x is the water velocity along the x-axis, 
v y is the water velocity along the y-axis, 
v z is the water velocity along the z-axis, 

V = *4 v x( a 4" a i) +v y( a 5' a 2) +v z( a 6* a 3)y c « tne expression for a correction to the phase corresponding to path 

differences of order v/c, where v is water velocity and c is the speed of sound in water, 

L is an upper limit to a summation, 

b n is a Legendre polynomial coefficient 

17 are ring geometry coefficients defined as follows: 



9 is the angle of measured energy away from the axis of wave propagation, 

m is taken to be zero, 

I is a dummy index being summed over, 

j is a dummy index being summed over, 

J is a Bessel function, 

P is a Legendre polynomial, and 

costi is the negative of the vertical leg away from the bistatic point. 

15. The system defined in Claim 9, wherein the closed form analytical expression for the spacial correlation function 
includes a Kirchhoff approximation for the scattering law function. 

16. The system defined in Claim 15, wherein the spatial correlation function is defined as follows: 



(l-m) I (j-m) ! r 
(1+m) ! (j+rn) I J 



sin^ cosd Pf(cos0) P?(cos6) dd 



where 



R = 



e 



x [ 



L 




(27r/r) 1/2 £ b x [P^cosp) J] 




where 



are ring geometry coefficients defined as follows: 




where 
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R 2 r 

o(0) = ^ duu J 0 (u) exp(-gu 2a ) 

4ir s±n 2 (6) cos 2 {8) J 



9 is the angle of measured energy away from the axis of wave propagation. 
R fo is a reflection coefficient, 
jo u is a variable over which the integral is taken, 

J is a Bessel function, 

Q = 2 (1. 2 a )/( ,2.2a) cos 2 (9)C 2 sin .2a (e) 

15 

Where 

q is an arbitrary variable chosen to simplify the equation 
k = 27t//u 

X is the wavelength of the propagated sound wave 
20 C h is a structure constant, also sometimes called a structure function parameter, and 

a is the effective relative angle of the azimuth. 

17. The system defined in Claim 9, wherein the closed form analytical expression for the spacial correlation function 
comprises a general function for band shapes. ...... 



18. The system defined in Claim 17, wherein the magnitude of the spatial correlation function is defined as follows: 



R peak 



2 .2 

30 {■: + (2 a -i)[2<x-x 0 yw] 2 } 1 ' a 



where 

35 R P eak is tne P eak spatial correlation, 

a is a wing parameter specifying the extent of the curve wing, 

X is a value on the curve, 

Xq is the center of the curve, and 

W is the width parameter of an angular scattering function defined as twice the distance from the center to the 
40 1/2 curve height point. 

19. A method for measuring velocities using a correlation sonar system having a transmitting transducer (208) and a 
plurality of receiving transducers (210). the method comprising the steps of transmitting a pulse of energy (110) 
in the direction of a reflective surface, and receiving the pulse energy echo (110') at each receiving transducer as 
45 a set of received values, the method further comprising maximum likelihood processing 1or estimating velocity 

using correlated values, 
characterised by: 

calculating correlation values from the received values associated vith multiple ones of the receiving trans- 
so ducers; 

providing initial estimates of a set of unknown parameters; and 

maximum likelihood processing the correlation values and initial estimates of parameters to maximize the 
probability of the system having a specific set of parameters given the correlation values. 

55 20. The method of Claim 1 9, additionally comprising the step of selecting a code for coding the pulse energy. 

21. The method of Claim 20, wherein the code is selected depending on relative measured velocity. 
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22. The method of Claim 20, wherein the code is selected depending on distance to the reflective surface. 

23. The method of Claim 19, wherein the energy is directed toward a plurality of reflective surfaces and the steps are 
repeated until a specific reflective surface is located. 

5 

24. The method of Claim 23, wherein the location of the specific reflective surface is by use of a matched filter. 

25. The method of Claim 19, wherein the set of unknown parameters includes a component of velocity. 

io 26. The method of Claim 19, wherein the maximum likelihood processing step includes a simplex algorithm. 

27. The method of Claim 1 9, wherein the maximum likelihood processing step includes evaluating the natural logarithm 
of the probability function at the initial estimate of parameters and at each parameter varied by an offset. 

is 28. The method of Claim 19, wherein the received values and correlation values are complex numbers. 
Patentanspruche 

20 1 . Korrelationssonarsystem (200) mit einem sendenden Wandler (209) zum Aussenden akustischer Energie auf eine 
reflektierende Oberflache zu, einer Impulserzeugungseinrichtung zum Erzeugen eines einzelnen Impulses (110) 
aus elektrischer Energie zu dem sendenden Wandler, einer Vielzahl von empfangenden Wandlern (210) zum Emp- 
fangen eines Impulsechos (110*) von der reflektierenden Oberflache, einer Vielzahl von Demodulatoren (222), 
wobei jeder Demodulator mit einem der empfangenden Wandler verbunden ist und einen Satz von Werten erzeugt, 

2S und einer MaximumwahrscheinKchkeitseinrichtung zum Schatzen der Geschwindigkeit unter Verwendung korre- 

lierter Werte, gekennzeichnetfernerdurcn das EinschlieRen einer Einrichtung (230) zum Autokorrelierender Werte 
fur jeden empfangenden Wandler und Kreuzkorrelieren zwischen verschiedenen empfangenden Wandlern mit 
einer Oder mehreren vorbestimmten Zeitverzogerungen, die kurzer als die Impulsdauer sind. 

30 2. System nach Anspruch 1 , wobei der Impuls eine Vielzahl codierter Elemente aufweist. 

3. System nach Anspruch 2, wobei jedes codierte Element eine Vielzahl von Tragerzyklen aufweist. 

4. System nach Anspruch 3, wobei jedes codierte Element rnehr als zwei und weniger als zehn Tragerzyklen aufweist. 

35 

5. System nach Anspruch 2, wobei die codierten Elemente phasencodiert sind. 

6. System nach Anspruch 5, wobei die codierten Elemente mit jeweils O-Grad- und 160-Grad-Phasen codiert sind. 

40 7. System nach Anspruch 2, wobei der Impuls Satze nichtidentischer Codeelemente umfaGt, so daR Nullen an jeder 
Seite einer vorgewahtten Autokorrelationsspitzenverzogerung vorgesehen werden. 

8. System nach Anspruch 1, wobei der Impuls eine zeitweilige Autokorrelationsspitze bei der vorbestimmten Zeit- 
verzogerung aufweist. 

45 

9. Korrelationssonarsystem (200) mit einer Einrichtung zum Erzeugen eines Impulses (110) in eine Flussigkeit mit 
reflektierenden Oberflachen, einer Vielzahl von Wandlern (210) zum Empfangen des Echosignals (110*), das durch 
den Impuls erzeugt wird, der auf die reflektierenden Oberflachen trifft, und mit einer Maximumwahrscheinlichkeits- 
einrichtung zum Maximieren der Wahrscheinlichkeit eines Satzes vorgewahlter Parameter, einschlieGlich zumin- 

50 dest einer Geschwindigkeitskomponente, mit Satzen von Korrelationen, die einen Parameterzustand zeigea wo- 

bei das System ferner gekennzeichnet ist durch das Aufweisen einer Korrelationseinrichtung zum Ausbilden von 
Korrelationen des Echosignals, das durch jeden Wandler mit einer vorbestimmten Zeitverzogerung empfangen 
wird, wobei das Echosignal durch den gleichen Wandler zu einer oder mehreren anderen Zeitverzogerungen emp- 
fangen wird und wobei die Echosignale durch die anderen Wandler mit einer Oder mehreren der anderen Zeitver- 

55 zogerungen emplangen werden, wobei die Korrelationen die raumliche Beziehung der Wandler bereitstellen. 

10. System nach Anspruch 9, wobei die Korrelationseinrichtung eine Einrichtung zum Berechnen von Korrelationen 
des Echosignals unter Verwendung einer raumlichen Korrelationsfunktion umfaGt, die definiert ist, wie folgt: 
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n = 0 



wobei 

b n ein Legendrepolynom-Koeffizient ist ( 

R„ eine raumliche Korrelationsfunktion ist, die den n-ten Term des Strahlmodells W(6) enthalt, 

W(9) das Systemstrahlmodell ist, und 

6 der Winkel der gemessenen Energie von der Achse der Wellenausbreitung weg ist. 

11. System nach Anspruch 10, wobei der Ausdruck fur R,, in der geschlossenen analytischen Form definiert ist, wie 
folgt: 



R B ° ^{l2l) x/2 n/2+n) i- J U2 . a (8) P n (cos€) 1 ^ 
6 ' a 2n<>l 



L<= F 



T n/2+k) i-* j l/a . k («) 



x P k (eose) (=1) 



(fl°l) ! i Jet 1 falls a 

n\ ! (Jc-l) t J [Jc(Jc+l) -n(a*l) J "«erade 1st 



» P„(eoo€) f-l ) ««>'*-n/» (Je-3.) 1 1 an falls a 

1st 



wobei 

n der n-te Term des Strahlmodells ist, 

\|/ = k[v x (a 4 -a 1 )+v y (a 5 -a 2 )+v z (a 6 -a 3 )]/c der Ausdruck fur eine Korrektur hinsichtlich der Phase ist, die zu Pha- 
sendifferenzen der GroBenordnung v/c gehort, wobei v die Wassergeschwindigkeit ist und c die Geschwin- 
digkeit des Schalls in Wasser ist, 

5 = (2n/X) x die gesamte Dreiachsenstrecke die Verzogerung vom bistatischen Punkt ist.. wobei X die Wellen r 
lange der sich ausbreitenden Schallwelle ist, 

J eine Besserfunktion ist, 

p ein Legendrepolynom ist, 

e ein Phasenwtnkel ist, der den relativen Beitrag einer vertikaten und horizontalen Versetzung beschreibt, 

und 

k = 2tz/X ist, wobei X die Wellenlange der sich ausbreitenden Schallwelle ist. 



12. System nach Anspruch 9, wobei der analytische Ausdruck in geschlossener Form fur die raumliche Korrelations- 
funktion eine Lambertsche Streugesetzfunktion umfaBt. 

13. System nach Anspruch 12, wobei die raumliche Korrelationsfunktion definiert ist, wie folgt: 
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R = _ f x [ 

2-0 

L m 

(2ir/r) 1 ' 2 £ b i tPj(cos^) £ i ' S J '°J J j.i/2< r > P^cosr,) * 
i • 
2 £ (cos/8) (-!)■ cos[jn(« 0 -a>] £ i * J Jf, v7^ 1/2 (r) P/fcosrj)]] 



wobei 

V = KVx{ a 4~ a i) +y/ yl a s~ a 2) + Vz( a G~ a 3)y c der Ausdruck fur eine Korrektur der Phase ist, die zu Streckenunterschie- 
den der GroBenordnung v/c gehort, wobei v die Wassergeschwindigkeit und c die Geschwindigkeit des 
Schalls im Wasser ist, 

L eine obere Grenze einer Summation ist, 

b n ein Legendrepolynom-Koeffizient ist, 

I™ Ringgeometriekoeffizienten sind, die definiert sind, wie folgt: 

IT 3 = (1/2+j) \Y m \ I i 3 '""! I f sind cosd p;(cosfl) P/Ccos*) d$, 
t-L+jn) : l j +m) i J 



wobei 

I ein Platzhafterindex ist, Ober den summiert wird, 

j ein Platzhafterindex ist, uber den summiert wird, 

m ein Platzhafterindex ist, uber den summiert wird, 
P ein Legendrepolynom ist, 

p der Nadirneigungswinkel des Meeresbodens ist, 

r = (2vJX) x der gesamten Dreiachsenstrecke die Verzogerung vom bistatischen Punkt (2tv x , 2tv y , 2tv 2 ) ist, 

wobei X die Wellenlange der sich ausbreitenden Schallwelle ist und t eine Zeitverschiebung ist, 
J eine Besselfunktion ist, 

cos ti das Negative der vertikalen Verzogerung von dem bistatischen Punkt weg ist, 

$ 0 -ot der Azimutwinkel ist : der zum Verhaltnis der wegfuhrenden Strecke hinsichtlich der y-Verzdgerung zur x- 
Verzogerung ist, die vom Azimutwinkel weg aufgenommen ist, unter dem der Meeresboden geneigt ist. 

14. System nach Anspruch 12, wobei die raumliche Korreiationsfunktion definiert ist, wie folgt: 



1-c 



wobei 

r = (2n/X) x die gesamte Dreiachsenstrecke der Verzogerung die Verzogerung vom bistatischen Punkt (2tv x , 
2tv y , 2tv z ) ist, wobei Xdie Wellenlange der sich ausbreitenden Schallwelle und t eine Zeitverschiebung ist, 

v x die Wassergeschwindigkeit langs der x-Achse : 
v y die Wassergeschwindigkeit langs der y-Achse. 
v z die Wassergeschwindigkeit langs der z-Achse ist, 

y = k(v x (a 4 -a 1 )+v y (a5-a2)+v z (a 6 -a 3 )]/c der Ausdruck fur eine Korrektur hinsichtlich der Phase zu Phasendiffe- 
renzen der Ordnung v/c gehorend ist, wobei v die Wassergeschwindigkeit und c die Geschwindigkeit des 
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Schalls im Wasser ist. 
L eine obere Grenze fur eine Summation ist, 

b n ein Legendrepofynom-Koeffizient ist, 

\™ Ringgeometriekoetfizienten sind, die definiert sind : wie folgt: 



10 

wobei 

9 • der Winkel der gemessenen Energie von der Achse der Wellenausbreitung weg ist, 

is 

m als Null angenommen wird, 

I ein Platzhalterindex ist, uber den summiert wird, 

j ein Platzhalterindex ist, uber den summiert wird, 

20 j eine Besselfunktion ist, 

P ein Legendrepoiynom ist, und 

costj das Negative der vertikalen Verzogerung vom bistatischen Punkt weg ist. 

15. System nach Anspruch 9, wobei der analytiscbe. Ausdruck in^geschlossener Form fur die raumliche Korrelations- 
25 f unktion eine Kirchhoffannaherung fur die Streugesetzfunktion umfaGt. 

16. System nach Anspruch 15, wobei die raumliche Korrelationsf unktion definiert ist, wie folgt: 

* • e - x c 

2]£ b, I? 0 P 1 (cos/3) 
(2Tr /^ )1/a E b > t*Mcos/3) £ HiJ^lr) PAcobti)* 
2 I Pi'(cos^) (-1)- cos [m (*,-<*>] £ *- J JT; ^. 1/a (r) P/Ccost,)]] 



40 wobei 



f™ Ringgeometriekoetfizienten sind : die definiert sind ; wie folgt: 



I? } = il./2+j) j j-jjjj j j j-g j jsind o <0) p;(cosfi) P/Ccos*) dd 



50 wobei 



R 7 

aid) » , _ '° — — duu J 0 {u) exp(-gu 2 «) 

4tt sin 2 (0) cosMd) 4 0 

9 der Winkel der gemessenen Energie von der Achse der Wellenausbreitung weg ist, 
R fo ein Reflexionskoeffizient ist, 
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u eine Variable ist, uber die integriert wird, 
J eine Bess elf unktion ist, 

q = 2 (1-2a, /t , 2 .2«, cos2(e)c 2 sin . 2a(e) 

wobei 

q eine befiebige Variable ist, die zum Vereinfachen der Gleichung gewahlt wurde, 
k = 2tl/X 

X die WeNenlange der stch ausbreitenden Schallwelle ist 

C h eine Strukturkonstante ist, die manchmai auch als ein Strukturfunktionsparameter bezeichnet wird, und 
a der effektive relative Winkel des Azimuts ist. 

17. System nach Anspruch 9, wobei der analytische Ausdruck geschlossener Form fur die raumliche Korrelations- 
funktion eine allgemeine Funktion fur Bandzustande ist. 

18. System nach Anspruch 17, wobei die Magnitude der raumlichen Korretationsf unktion definiert ist, wie folgt: 

a "Spftze 

n- 2 2 

{1 + (2 a -1)[2(X-X 0 )/VV] 2 } 1/a 

wobei 

^spitze d ' e raumliche Spitzenkorrelation ist, 

a ein Schenkelparameter ist : der das AusmaG des Kurvenschenkels spezifiziert, 

x ein Wert der Krummung ist, 

Xq die Mitte der Krummung ist, und 

W der Breitenparameter einer Winkelstreuf unktion ist, der als das Zweifache der Strecke von der Mitte 

zum 1/2-Kurvenhohenpunkt definiert ist. 

19. Verfahren zum Messen von Geschwindigkeiten unter Verwendung eines Korrelationssonarsystems mit einem sen- 
denden Wandler (208) und einer Vlelzaht von empfangenden Wandlern (210), wobei das Verfahren die Schritte 
des Sendens eines Energieimpulses (1 1 0) in der Richtung einer reflektierenden Oberflache, und des Empfangens 
des Impulsenergieechos (110') bei jedem empfangenden Wandler als einen Satz empfangener Werte aufweist, 
wobei das Verfahren ferner eine Maximum-Likelihood-Verarbeitung zum Abschatzen der Geschwindigkeit unter 
Verwendung korrelierter Werte, 

aufweist, 

gekennzeichnet durch 

das Berechnen von Korrelationsverten aus den empfangenen Werten, die zu vielen der empfangenden Wand- 
ler gehoren; 

Bereitstellen anfanglicher Abschatzungen eines Satzes unbekannter Parameter; und 
eine Maximum-Likelihood- Verarbeitung der Korrelationswerte und der anfanglichen Abschatzungen von Pa- 
rametern zum Maximieren der Wahrscheinlichkeit des Systems, das einen bestimmten Satz von Parametern 
aufweist, die den Korrelationswerten zugeordnet sind. 

20. Verfahren nach Anspruch 1 9, zusatzlich den Schritt des Auswahlens eines Codes zum Codieren der Impulsenergie 
aufweisend. 

21. Verfahren nach Anspruch 20, wobei der Code abhangig von der relativ gemessenen Geschwindigkeit ausgewahlt 
wird. 

22. Verfahren nach Anspruch 20, wobei der Code abhangig von der Strecke zur reflektierenden Oberflache ausgewahlt 
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wird. 

23. Verfahren nach Anspruch 19, wobei die Energie zu einer Vielzahl reflektie render Oberflachen hin gerichtet wird 
und die Schritte wiederholt werden, bis eine bestimmte reflektierende Oberflache lokalisiert ist. 

24. Verfahren nach Anspruch 23, wobei die Lage der bestimmten reflektierenden Oberflache durch Anwenden eines 
angepaBten Fitters bestimmt wird. 

25. Verfahren nach Anspruch 1 9, wobei der Satz unbekannter Parameter eine Geschwindigkeitskomponente umfaBt. 

26. Verfahren nach Anspruch 19, wobei der M^imum-Likelihood-Verarbeitungsschritt einen Simplexalgorithmus um- 
faBt. 

27. Verfahren nach Anspruch 1 9, wobei der Maximum-Likelihood-Verarbeitungsschritt das Berechnen des naturlichen 
Logarithmus der Wahrscheinlichkeitsfunktion bei der anfangiichen Abschatzung von Parametem und bei jedem 
Parameter durch einen Offset variiert, aufweist. 

28. Verfahren nach Anspruch 19, wobei die empfangenen Werte und Korrelationswerte komplexe Zahlen sind. 
Revendications 

1. Systeme sonar a correlation (200), qui possede un transducteur d'emission (208) servant a'emettre une energie 
acoustique sur une surface reflechissante,..unjnoyen generateur d'impulsipn servant a produire. une unique im- 
pression (110) d'energie electrique a destination du transducteur d'emission, une pluralite de transducteurs de 
reception (210) servani a recevoir un echo d'impulsion (110*) de la part de la surface reflechtssante, une pluralite 
de demoduiateurs (222), chaque demodulateur etant connecte a Tun des transducteurs de reception et produisant 
un ensemble de valeurs, et un moyen d'application de probability maximale permettant d'estimer la Vitesse a I'aide 
de valeurs correlees, caracterise en ce qu'il comporte en outre : 

un moyen (230) servant a determiner par autocorrelation tesdttes valeurs pour chaque transducteur de re- 
ception, et a effectuer une intercorrelation entre differents transducteurs de reception, a un ou plusieurs temps de 
retard predetermines moindres que la largeur d'impulsion. 

2. Systeme selon la revendication 1, ou I'impulsion comprend une pluralite d'elements codes. 

3. Systeme selon la revendication 2 f ou chaque element code comprend une pluralite de cycles de porteuses. 

4. Systeme selon la revendication 3, ou chaque element code comprend plus de deux et moins de dix cycles de 
porteuses. 

5. Systeme selon la revendication 2, ou les elements codes sont codes en phase. 

6. Systeme selon la revendication 5, ou les elements codes sont codes avec des phases de 0° ou bien de 180°. 

7. Systeme selon la revendication 2, ou Timpulsion comporte des ensembles d'elements de code non identiques de 
facon a produire des zeros de part et d'autre d'un retard de crete d'autocorrelation preselectionne. 

8. Systeme selon la revendication 1, ou I'impulsion possede une crete d'autocorrelation temporelle au temps de 
retard predetermine. 

9. Systeme sonar a correlation (200), qui possede un moyen servant a produire une impulsion (110) dans un liquide 
possedant des surfaces reflechissantes, une pluralite de transducteurs (210) servant a recevoir le signal d'echo 
(110') cree par I'impulsion venant en contact avec les surfaces reflechissantes, et un moyen d'application de pro- 
bability maximale servant a maximiser la probability d'un ensemble de parametres preselectionnes, comportant 
au moins une composante de Vitesse, ayant des ensembles de correlations indicatifs d'un etat de parametres, 
ledit systeme etant caracterise en ce qu'il comprend en outre : 

un moyen de correlation servant a former des correlations du signal d'echo recu par chaque transducteur a 
un temps de retard predetermine avec le signal d'echo recu par le meme transducteur a un ou plusieurs autres 
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temps de retard, et avec les signaux cf echo recus par les autres transducteurs k un ou plusieurs des autres temps 
de retard, lesdites correlations permettant la mise en relation spatiale des transducteurs. 

10. Systeme selon la revendication 9, ou le moyen de correlation comporte un moyen permettant de calculer des 
s correlations du signal d'echo a Paide d'une fonctton de correlation spatiale dSfinie comme suit : 



10 
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40 
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SO 



n = 0 



ou: 

b n est un coefficient de polynome de Legendre, 

R n est une lonction de correlation spatiale faisant intervener le n feme terme du diagramme de faisceau W (8), 
W(9) est le diagramme de faisceau du systeme, et 

6 est Tangle d'£nergie mesuree par rapport a I'axe de propagation des ondes. 
11. Systeme selon la revendication 10, ou I'expression analytique de forme fermee pour est definie comme suit : 



25 R n = c'V {[^~ Y* (l/2+»)i- W5) />/,(cose) 



1 



5 J ' ' v ~ ' 2/;+l 

\ J ~Y a T] (1/2+*) i^iWS) 

5 k= 1 

impair 

(//-!)!!*!! 

x P,. (cos £)(-!) V2 si nest pair 

/;!! !! [k(k+] )-//(//+!)] 



oo 

5 k = 0 

pair 

(k-l)WnU 

x P t (cose) (-1) ("* fc - lv 2 si n est impair } 

k [n(n+\yk(k+\y\ 



ou: 



n est le n i6me terme du diagramme de faisceau, 

\|/ = k [v x (a 4 - a, ) + v y (a 5 - a 2 ) + v 2 (a 6 - a 3 )] /c, qui est I'expression relative a la correction apport6e a la phase 
correspondant a des differences de trajet d'ordre v/c t ou v est la Vitesse de I'eau et c est la Vitesse du son 
55 dans Feau, 

8 = (2idX) x la distance totale, suivant les trois axes, pour le retard par rapport au point bistatique ( ou X est la 
longueur d'onde de I'onde sonore propagee, 
J est une fonctton de Bessel, 



27 



EP 0 619 024 B1 



P est un polynome de Legendre, 

e est un angle de phase decrivant la contribution relative du deplacement vertical et horizontal, et 
k = 2iUk, ou X est la longueur d'onde de Tonde sonore propagee. 

1 2. Systeme selon la revendication 9, ou Texpression analytique sous forme fermee relative a la fonction de correlation 
spatiale comporte une fonction de loi de diffusion de Lambert. 

13. Systeme selon la revendication 12, ou la fonction de correlation spatiale est definie comme suit : 

/?= X 

L 

2£ MfoP^cosp) 
1 = 0 

L oo 
[<2n/r)'« 2 h \ ^i( C0S P) Z [i r ?J W < r > P j< cos7 l> + 

1=0 j=0 



1 00 

2^ P; (cos(J) •(.!)« cos [w( 0 o -a)] £ i^7/ W W ~VJ (cost,)]] 
w ~ 1 7 = m 

ou: 

\j/ = k [v x (a 4 - a, ) + v y (a 5 - a 2 ) + v z (a 6 - a 3 )]/c, qui est Texpression relative a une correction apportee a la 
phase correspondant a des differences de trajet d'ordre v/c, ou v est (a Vitesse de I'eau et c est la Vitesse du son 
dans I'eau, 

L est une limite superieure pour la sommation, 
b n est un coefficient de polynome de Legendre, 

I™ designe des coefficients de geometrie annulaire definis comme suit : 

=(i/2t/) ! °" m) ! [sine cose />7(cose)p m (cose) de - 

e, 



I est un indice fictif sur lequel la semination est effectuee, 
j est un indice fictif sur lequel la sommation est effectuee, 
m est un indice fictif sur lequel la sommation est effectuee, 



P est un polynome de Legendre, 

P est ('angle d'inclinaison du fond marin par rapport au nadir, 

r = (2-nJX) x la distance totale, par rapport aux trois axes, pour le retard par rapport au point bistatique (2tv x , 
2tv y , 2tv z , ou X est la longueur d'onde de I'onde sonore propagee, et t est un decalage tempore!, 
J est une fonction de Bessel, 

cosri est la vaieur negative du retard vertical par rapport au point bistatique, 

0 o -a est Tangle azimutal correspondant au rapport de la distance d'ecartement sur le retard y a la distance 
d'ecartement sur le retard x, indexe par rapport a Tangle azimutal suivant lequel le fond de la mer est incline. 
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14. Systeme selon la revendication 12, ou la lonction de correlation spatiale est de7inie comme suit : 



(2^/r) I ^e^ / 



L 



OC 



R = 



L 



x [2 



*i X HI?; i WW Pj(cosn)] 



2 s a./; 



ou : 



r = (2tc/X) x la distance totale, suivant les trois axes, pour le retard par rapport au point bistatique (2tv x , 2xv y , 
2tv z ), ou X est la longueur d'onde de I'onde sonore propagee : et x est un retard temporel, 



\y = \|/ = k [v x (a 4 - a,) + v y (a 5 - a 2 ) + v z (a 6 - a 3 ))/c, qui est 1'expression relative a la correction apportee a la 
phase correspondant a des differences de trajet d'ordre v/c, ou v est la Vitesse de I'eau et c est la Vitesse 
du son dans I'eau, 



9 est Tangle de I'energie mesuree par rapport a I'axe de la propagation de I'onde, 
m est pris egal a zero, 

I est un indice fictif sur lequel la sommation est effectuee, 
j est un indice fictif sur lequel la sommation est effectuee, 

J est une fonction de Bessel, 

P est un polynome de Legendre, et 

cost] est la valeur negative du retard vertical par rapport au point bistatique. 

15. Systeme selon la revendication 9, ou 1'expression analytique sous forme fermee relative a la fonction de correlation 
spatiale comporte une approximation de Kirchhoff pour la fonction de loi de diffusion. 

16. Systeme selon la revendication 15, ou la fonction de correlation spatiale est definie comme suit : 



v x est la Vitesse de I'eau suivant I'axe x, 
v y est la Vitesse de I'eau suivant I'axe y, 
v z est la vitesse de I'eau suivant I'axe z, 



L est une limrte superieure pour la sommation, 

b n est un coefficient de polynome de Legendre, 

I? sont des coefficients de geometrie annulaire definis comme suit : 




e, 



ou 
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L 

2£ MfoP^cosP) 
1 = 0 



L oo 

[(2K/r) 1/2 X) *i ^i( C0S P)Z H 1 i j h+i* W P i (cosTl) + 
1=0 j= 0 



1 oo 
2^ P 7 < C0S P> cos [ifi( 0 o -a)] X HI"! WW ^ (wsn)]] 



ou 

I™ sont des coefficients de geometrie annulaire definis comme suit : 

I" = (1/2 +y) ' Q ' m) ' f sin B o (6) P (cos9 ) /> " (cosG ) d9 

(I+w) ! (j+w) ! J J 



ou 



, oo 



r,: f 
0(9)= — — rrr; — r; J du u J o < u) ex P H" 2 * > 

4n sin 2 (9) cos 2 (9) J 



9 est Tangle de Tenergie mesuree par rapport a I'axe de propagation de I'onde, 
R fo est un coefficient de reflexion, 

u est une variable sur laquelle I'integration est effectuee, 
J est une fonction de BesseL 

q = 2( v2 °) k( 2 ' 2a ) cos 2 (0) C 2 sin' 2a (0) 



q est une variable arbitrage choisie pour simplifier Pequation, 

X est la longueur d'onde de I'onde sonore propagee, 

C h est une constante de structure, queiquefois appelee parametre de fonction de structure, et 
a est Tangle relatif effectif de I'azimut. 

17. Systeme selon la revendication 9, ou ('expression analytique sous forme fermee pour la fonction de correlation 
spatiale comprend une fonction generale pour des formes en bandes. 

18. Systeme selon la revendication 17, ou ('amplitude de la fonction de correlation spatiale est definie comme suit : 
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est la correlation spatiale de crete, 

est un parametre de section specifiant I'etendue de la section de courbe, 
est une valeur sur la courbe, 
est le centre de la courbe, et 

est le parametre de largeur d'une fonction de diffusion angulaire definie comme deux fois la distance du 
centre au point a ni-hauteur de la courbe. 

19. Procede de mesure de vitesses a I'aide cfun systeme sonar a correlation possedant un transducteur d'emission 
(208) et une pluralite de transducteurs de reception (210), le procede comprenant les operations consistant a 
emettre une impulsion d'energie (1 1 0) dans la direction d'une surface reflechissante, et a recevoir Techo d'energie 
d'impulsion (110*) en chacun des transducteurs de reception sous la forme d'un ensemble de valeurs recues, le 
procede comprenant en outre un traitement d'application de probabilite maximale permettant d'estimer la vitesse 
a I'aide de valeurs correlees, 

caracterise par les operations suivantes : 

calculer des valeurs de correlation a partir des valeurs recues associees a plusieurs des transducteurs de 
reception ; 

produire des estimations initiales d'un ensemble de parametres inconnus ; et 

effectuer un traitement d'application de probabilite maximale cur les valeurs de correlation et des estimations 
initiales de parametres afin de maximiser la probabilite du systeme d'avoir un ensemble particulier de para- 
metres etant donne les valeurs de correlation. 

20. Procede selon la revendication 19, comprenant en outre ('operation qui consiste a selectionner un code afin de 
coder I'energie d'impulsion. 

21. Procede selon la revendication 20, ou on choisit le code en fonction de la vitesse mesuree relative. 

22. Procede selon la revendication 20, ou on choisit le code en fonction de la distance a la surface reflechissante. 

23. Procede selon la revendication 19 t ou I'energie est dirigee vers une pluralite de surfaces reflechissantes et on 
repete les operations jusqu'a ce qu'une surface reflechissante particuliere soil localisee. 

24. Procede selon la revendication 23, ou la localisation de la surface reflechissante particuliere est faite par utilisation 
d'un filtre apparie. 

25. Procede selon la revendication 19, ou I'ensemble de parametres inconnus comporte une composante de vitesse. 

26. Procede selon la revendication 19, ou I'operation d'application de probabilite maximale comporte un algorithme 
simplex. 

27. Procede selon la revendication 19, ou I'operation d'application de probabilite maximale comporte revaluation du 
logarithme naturel de la fonction de probabilite pour I'estimation inrtiale des parametres et pour chaque parametre 
ayant subi un decalage. 

28. Procede selon la revendication 1 9, ou les valeurs recues et les valeurs de correlation sont des nombres complexes. 
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